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We compute the fractal dimension of clusters of inertial particles in mixing flows at finite values 
of Kubo (Ku) and Stokes (St) numbers, by a new series expansion in Ku. At small St, the theory 
includes clustering by Maxey's non-ergodic 'centrifuge' effect. In the limit of St ^ oo and Ku 
so that Ku^St remains finite) it explains clustering in terms of ergodic 'multiplicative amplification'. 
In this limit, the theory is consistent with the asymptotic perturbation series in [Duncan et a/., Phys. 
Rev. Lett. 95 (2005) 240602]. The new theory allows to analyse how the two clustering mechanisms 
compete at finite values of St and Ku. For particles suspended in two-dimensional random Gaussian 
incompressible flows, the theory yields excellent results for Ku < 0.2 for arbitrary values of St; the 
ergodic mechanism is found to contribute significantly unless St is very small. For higher values 
of Ku the new series is likely to require resummation. But numerical simulations show that for 
Ku ~ St ~ 1 too, ergodic 'multiplicative amplification' makes a substantial contribution to the 
observed clustering. 
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Introduction. The dynamics of independent particles 
in complex mixing flows is a problem of fundamental im- 
portance in the natural sciences, and in technology. The 
motion of the particles is commonly approximated by 

r = 7[u(r, t) — v] . (1) 

Here r is the position of a suspended particle, and v = r 
is its velocity. Dots denote time derivatives, 7 is the 
rate at which the inertial motion is damped relative to 
the fluid, and u{r^t) is the velocity of a randomly mix- 
ing or turbulent incompressible flow. It is a surprising 
fact that even though u is incompressible, the suspended 
particles may nevertheless cluster together [T. The ef- 
fect is illustrated in Fig. [l^,b. Possible consequences 
of this phenomenon have been discussed in a wide range 
of contexts: rain initiation from turbulent clouds |2j-t4 , 
grain dynamics in circumstellar accretion disks 0|6], and 
plankton dynamics [7 , to name but a few. 

Despite its significance, clustering of particles in mix- 
ing flows is still not well understood. Two very different 
explanations of the phenomenon have been put forward. 
Maxey fP^ discussed the problem in the limit of small in- 
ertia, corresponding to small values of the 'Stokes num- 
ber' St = (7r)~^. Here r is the relevant characteristic 
time scale of the flow (the Kolmogorov time in turbulent 
flows, for example). For < St <C 1, the particles are 
argued to be 'centrifuged' out of regions of high vorticity 
of u{r^t). The approach rests on instantaneous correla- 
tions between particle positions and fluid velocities. It 
has been refined by many authors jSHTO] and is com- 
monly used to interpret results of experiments [HI [12] , 
and of direct numerical simulations [13]. But the 'cen- 
trifuge' mechanism relies on a small-St expansion, while 
clustering in turbulent flows is observed to be strongest 
and thus of most interest at St ~ 1. 

A very different clustering mechanism was proposed 
[T4] in the limit of large St and small 'Kubo numbers'. 




FIG. 1: a Clustering of particles in a two-dimensional 
incompressible flow it(r,t) = V A '0(r,t)e3. Here 63 is 
the unit vector ± to the plane. The Gaussian random 
function '0(r,t) satisfies (ip) = and ('0(r, t)'0(0, 0)) = 
(no^^/2) exp [— |r|^/(2?7^) — |t|/T]. Green contours corre- 
spond to high vorticity of it(r,t), blue contours to high 
strains. Particle number density: white (low density) to red 
(high density). Parameters: Ku = 0.1, St = 10, t = 165t. b, 
same but for Ku = 10, St = 0.025, t = 0.32r. c Param- 
eter plane for inertial particles in mixing flows. Region 1: 
clustering is caused by ergodic 'multiplicative amplification', 
see text. Region 2: the non-ergodic 'centrifuge' mechanism 
is important. Turbulent flows correspond to Ku ^ 1, strong 
clustering is observed for St ^ 1, region 3. In region 4 the 
new perturbation expansion is accurate (schematic). 

The Kubo number [Mj [15] Ku = uqt/t] characterises 
fluctuations of u{r,t) {uq and r] being its characteristic 
velocity and length scales). In the limit St ^ 00 (and 
Ku ^ so that = Ku^St/2 remains constant), the 
particles experience the velocity field as a white-noise sig- 
nal, and sample it in an ergodic fashion: the fluctuations 
of u(rt, t) (and its derivatives) along a particle trajectory 
Tt are indistinguishable from the fluctuations of u{ro,t) 
at the fixed position tq. This case corresponds to region 
1 in the phase diagram Fig. [ij^, and in this limit the in- 
stantaneous configuration of u{rt^t) is irrelevant to the 
dynamics of the suspended particles. But they may nev- 
ertheless cluster by the mechanism of 'multiplicative am- 
plification': small line-, area-, and volume elements ran- 
domly expand and contract. Depending upon whether 
the random product of expansion and contraction fac- 
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tors increases or decreases as t oo, one may observe 
fractal clustering in region 1. The fractal dimension di, 
is determined by the history of these factors. It can be 
computed in terms of 'Lyapunov exponents' [T4( [T6]. 

Figs. [l^,b show both mechanisms at work: at large 
values of St, (region 1 in Fig. [l]^) there is no discernible 
influence of the instantaneous u{rt^t) on the particle dis- 
tribution. At small St, (region 2 in Fig. [T]::), by contrast, 
the particles are seen to avoid regions of high vorticity 
(Fig. [T}3, similar to Fig. 8 in ^j). In short, in limit- 
ing cases (regions 1 and 2 of Fig. [T]^) the mechanisms 
of clustering are understood. But how ergodic and non- 
ergodic effects compete in the major part of the phase 
diagram Fig. is not known (in particular not for the 
experimentally most relevant region 3 where Ku and St 
are of order unity). In [TF non-ergodic effects were char- 
acterised by correlating the degree of clustering with the 
probability of particles avoiding rotational regions of the 
flow. The interpretation of these numerical results, how- 
ever, is complicated by the fact that this probability is 
significantly enhanced even when clustering is weak. In 
order to understand the importance of non-ergodic and 
ergodic effects, an analytical theory is required, valid at 
finite Stokes and Kubo numbers. 

Summary. Here we derive a perturbation expansion 
for the Lyapunov exponents of particles in random flows, 
valid at finite St and Ku. We compute the Lyapunov 
fractal dimension c^l, and characterise fractal clustering 
in terms of the 'dimension deficit' Al = d — di,. For par- 
ticles suspended in two-dimensional random Gaussian in- 



compressible flows, the new theory yields reliable results 
for Ku < 0.2 for arbitrary values of St, region 4 in Fig. [T]^ 
(schematic). We find, first, that for small values of St, the 
'centrifuge' mechanism dominates, and Al = 6Ku^St^, 
consistent with [20 , 21^. Second, in the limit of Ku ^ at 
finite values of St, non-ergodic effects remain important. 
Third, in region 1 of Fig. clustering is found to be 
entirely due to ergodic 'multiplicative amplification' [14], 
and Al = 12e^ cx St. Fourth, in general we find that the 
ergodic mechanism contributes substantially to cluster- 
ing, unless St is very small. Fifth we show by numerical 
simulations of the model that at Ku ~ 1 and St ~ 1, er- 
godic 'multiplicative amplification' makes a substantial 
contribution to the observed clustering. 

Method and results. Eq. ([T]) cannot be explicitly solved, 
since u depends upon the particle position at time t. 
The implicit solution of Eq. ([T]) becomes (dimensionless 
variables r' = r/ry, t' = t/r, v' = v/uq, and u' = u/uq) 



Srt = rt-ro = Ku 
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Here and in the following the primes are omitted. We 
seek an approximate solution by expanding u{rt^t) in 
powers of Srt. Since according to Eq. Srt is of order 
Ku, iteration generates an expansion of n(rt, t) in powers 
of Ku. To second order, for example, we find 
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Here Greek indices denote the components of n, and for 
our purposes can be set to zero. The coefficients in 
(|3| are expressed in terms of u and its derivatives at 
the fixed position ro, with known statistical properties. 
This procedure can in principle be extended to any or- 
der in Ku, but in practice it is limited by the number 
of nested integrals appearing in (|3| for higher orders. A 
C-program was written to symbolically evaluate the inte- 
grals. One may expand other functionals of the particle 
trajectories, such as the strain matrix A(rt,t) with ele- 
ments Aq,^ = duo^/dvj^. Previous analytical results on 
the clustering of inertial particles [Ml [19] rest on the 



'ergodic assumption' that the distribution of the strain 
matrix A(r^, t) at the particle position rt can be approxi- 
mated by its distribution at ro- This is satisfied in region 
1, but what are the corrections outside this region? For 
the two-dimensional incompressible (TrA = 0) random 
flow described in Fig. [l] we find: 



TrA2 = / lim -/ dtTrA'(n,t) 

XT^ooTJo 

2Ku'^St(4 + 52St + 293St' 



6Ku'St 



(1- 
548St^ - 



-St)2(l + 2St) 
- 297St'^) 



(1 + St)4(2 + St)(l + 2St)2(l + 3St) 



(4) 



3 



to order Ku^. The average in Q consists of a long-time 
average along the particle trajectory r^, and an average 
over initial conditions tq (denoted by (•••)). At finite 
values of the Kubo number, TrA^ differs from its ergodic 
average (which vanishes in incompressible fiows): the dy- 
namics is not strictly ergodic. In compressible fiows, we 
find TrA ^ (the ergodic average still vanishes). This 
is consistent with a result [22^ for the average strain in 
the advective limit of a one-dimensional (compressible) 
model. It was shown in [22 that the average strain must 
be taken into account to obtain the known result [9 for 
the advective Lyapunov exponent in this model. 

The question is now: how does non-ergodicity affect 
the spatial distribution of the particles? The latter is 
characterised by the Lyapunov exponents of the parti- 
cle fiow, obtained by linearising Eq. ([T]). The maximal 
exponent Ai is given by 

1 

Ai = Ku hm - / dtni{rt,t) ■Z{rt,t)ni{rt,t) , (5) 
Z = St"^(A-Z)-KuZ^ ni=Ku(n2 •Zni)n2. (6) 
Here rii {712) is the unit vector in the Sr ((5r)-direction, 



and Z is the matrix with elements Za/3 = dvc^/dr^. In 
region 1, Eqs. \b\o\ were solved in [14[[T9]. At finite val- 
ues of Ku, we compute Ai by generalising the procedure 
that led to Eq. ([s]). Starting from the implicit solution 
([2| of ([1]), and the implicit solutions of (|6|: 

Z{rt,t) = e-*/s*Z(ro,0)+ / dtie"^ [A(n, , ti)/St 

-KnZ{rt,Mf], (7) 

ni{rt,t) = ni(ro,0) +Ku / dti[n2{rt^,ti) ■ Z{rt^,ti) 

Jo 

xni(rt,,ti)]n2(n,,ti) , (8) 

we expand A(r^,t), Z{rt,t), ni(rt,t), and n2{rt^t) in 
powers of drt. Iterating and averaging along particle 
trajectories as well as over initial conditions yields an 
expansion of Ai in powers of Ku, with St-dependent co- 
efficients. The sum Ai + A2 = KuTrZ is computed in 
a similar fashion. For particles in a two-dimensional in- 
compressible random Gaussian ffow (c.f. Fig. fl]) we find 
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to order Ku'. This is our main result. Eq. ^ yields the 
ergodic expansion [1£ Xi/{jr) = 2e^ - 206^^ + 4806^ + . . . 
in region 1. 



-GKu'^St^ + 4Ku'St^ + 0(Ku^ 



As St 0, Eq. (10) reflects Maxey's non-ergodic cen- 
trifuge mechanism: According to ([T]), a particle is ad- 
vected by an effective velocity field v with compressibil- 
ity V • i; = TrZ. Maxey's result [l] |2] is obtained by 
expanding Z ^ Z^^) + Z^^^St in Eq. One finds: 

V • v = -KuSt TrA^lg^^Q (note that A(rt,t) depends 
upon the Stokes number because the particle trajectory 
rt depends upon St). This result shows that particles 
tend to aggregate in regions of high strain or low vor- 
ticity. However, since the velocity field is homogeneous, 
this lowest-order term averages to zero in incompress- 
ible flows. Expanding Z to second order in St, one finds 

Inserting Eq. ( 4 ) yields 



V'V = -KuSt^ SstTr(A2) 



St=o 



limit of (10). 



Ai + A2 = Ku V • i; = 
consistent with the St 

We conclude that our new expansion ( |9|10D of the Lya- 
punov exponents correctly describes the different clus- 
tering mechanisms in the advective and ergodic regions 
(Fig. [1]::): the 'centrifuge' effect and ergodic 'multiplica- 
tive amplification'. More importantly, Eqs. ( 9|10 ) al- 
low to determine how the relative importance of the two 
mechanisms depends on Ku and St, as follows. 

In two-dimensional incompressible fiows, the fractal di- 
mension deficit is given by Al = (Ai + A2)/A2. Fig. [2^ 
shows that the new theory explains the limiting cases in 
region 1 in the parameter plane (Al = 6 Ku^ St = 12 6^), 
and in region 2 (Al = 6Ku^St^). The theory also ex- 
plains the cross-over between these two behaviours and 
compares well with results of numerical simulations. In 
order to further increase the accuracy, more terms than 
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FIG. 2: a Fractal dimension deficit as a function of e = 
Ku^St/2. Numerical simulations of the model described in 
Fig. [l](Ku = 0.02 (0), Ku = 0.05 (A ) Ku = 0.1 (o), and 
Ku = 0.2 (v); theory according to Eqs. ( 9|10 ), solid lines, and 
limiting behaviours Al oc St and Al oc St (dashed lines), 
b Same but for larger values of = Ku^St/2. Numerics, 



Ku = 0.1 (o) and Ku = 1 (□); theory, Eqs. ( |9jl0| ), Ku = 0.1 
(solid line), c Relative importance of non-ergodic and ergodic 
contributions (see text), symbols and parameters as in b. 



computed in Eqs. ( 9|10 ) must be included. The se- 
ries is likely to be asymptotic requiring re-summation, 
and there may be additional non-analytic contributions 
[HI [19] . We note that for a slightly different estimate of 
the fractal dimension (the 'correlation dimension'), the 
lowest-order behaviour of the dimension deficit in region 
2, Ac oc St^, was computed using different methods in 
P fTOl 1231425;. For the model described in Fig.[lJ we find 
Ac = 12 Ku^St^ to lowest order in Ku and in St in region 
2, and Ac = 24 in region 1 (see also [26j). 

These results raise the question: how important are 
non-ergodic contributions to fractal clustering at larger 
values of St and at finite Kubo numbers? The answer is 
summarised in Fig. [2j3,c showing the dimension deficit 
Al compared to an ergodic approximation, A^^^. The 
latter incorporates finite-time correlations of the velocity 
field u but neglects non-ergodic effects. Ergodic approx- 
imations to the Lyapunov exponents and to Al (referred 
to as A^^l ^L^) obtained by expanding Eqs. (5 6) 
as before, but replacing A(rt, t) in ([6| with A(ro, t). Ihe 
resulting analytical expressions for A^^f determined 
by the fluctuations of A(ro, t). This is in contrast to Ai^2, 
Eqs. (|9|10[), which are determined by the fluctuations of 



n(ro,t) and its derivatives. The ergodic approximation 
allows for finite values of Ku and St but must fail in 
the limit St 0, since the 'centrifuge' mechanism is not 
accounted for. In particular, to lowest order in Ku the 
exponents A^^^ and A2^^ are found to depend upon St. 
The exact exponents ( 9|10 ) by contrast, are independent 
of St to lowest order in Ku: Ai = Ku^ and A2 = — Ku^. 
This implies in particular that earlier results for the Lya- 
punov exponents in region 1 [IT, T9^ are in fact exact to 
lowest order in Ku for arbitrary values of St. This is due 
to the cancellation of two errors: neglecting non-ergodic 
effects, and neglecting finite-time correlations. 

Non-ergodic effects dominate the clustering when 
(Al — Al^^)/Al is close to unity (they are negligible 
when this ratio is close to zero). We have determined 
A^^l ^ order Ku^ [27^ . We have also performed 



computer simulations of this 'ergodic model' by calculat- 
ing monodromy matrices (Eq. (30) in [21]) with A(ro,t) 
evaluated at the fixed position vero. Fig.^ shows that 
for Ku = 0.1, non-ergodic effects dominate at small val- 
ues of but are negligible when clustering is largest, near 
the peak in Al shown in Fig. [2)3. Also shown are the an- 
alytical result for Al derived from Eqs. ( |9|10 ), and the 
corresponding expansion of (Al — Al^^)/Al. We observe 
good agreement. For Ku = 1, by contrast, the first terms 
in the perturbation expansion ( |9|10[ ) do not give reliable 
results. But computer simulations show that non-ergodic 
effects are present for the whole range of Stokes numbers 
displayed in Fig. However, Fig. ^ also clearly shows 
that both mechanisms contribute in region 3. We em- 
phasise that ergodic clustering by 'multiplicative ampli- 
fication' makes a substantial contribution in this region, 
(Al-A^l'^)/Al^0.3. 
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